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We apply the maximum entropy method to extract the spectral functions for pseu- 
doscalar and vector mesons from hadron correlators previously calculated at four dif- 
ferent lattice spacings in quenched QCD with the Wilson quark action. We determine 
masses and decay constants for the ground and excited states of the pseudoscalar and 
vector channels from position and area of peaks in the spectral functions. We obtain 
the results, m^^ = 660(590) MeV and rup^ = 1540(570) MeV for the masses of the 
first excited state masses, in the continuum limit of quenched QCD. We also find un- 
physical states which have infinite mass in the continuum limit, and argue that they 
are bound states of two doublers of the Wilson quark action. If the interpretation is 
correct, this is the first time that the state of doublers is identified in lattice QCD 
numerical simulations. 

11.15.Ha, 12.38.Gc 



I. INTRODUCTION 



The spectral function of hadron correlation functions contains information not only on the mass 
of the ground state but also other quantities such as the masses for excited states, and decays and 
scatterings of hadrons. In lattice QCD simulations one can numerically obtain an euclidean time 
correlation function D{t) of an operator 0(r), which is related to the spectral function /(cj) of this 
correlator through 

i?(r) = (0|O(r)Ot(0)|0) 

= j dujf{L0)K{L0,T), (1) 

where K{t,u}) is a kernel of the Laplace transformation given by 

for < r < T with the periodic boundary condition, where T is the lattice size in the euclidean time 
direction. A typical form of /(w) is 

f[u)^ Za5[u^Eo) + J[u-u>2mo), (2) 

where Eq is the energy of the ground state \Eq) coupled to the operator O and Zq = |(0|0|i?o)p, and 
/(cj) represents the continuous spectrum which starts at w = 2mo for the 2-particle state. 

In principle one can extract all the information for the states which can couple to the operator O 
from the spectral function ,f{uj). In a usual analysis of lattice QCD simulations, however, only the 
mass (or energy) of the ground state _Bo and its amplitude can be reliably extracted from the 
asymptotic behavior of the correlation functions at large euclidean times, 

D{t) Zoe-^'"', r ^ oo. 



1 



Numerically it is unstable to extract masses of excited states with a multi-exponential fit, so that a 
simultaneous fit to several correlation functions which have the same intermediate states with different 
amplitudes becomes necessary to stabilize the result. Similar but more difficult problems appear in 
the calculation of the decay amplitude [Qj^ . 

If one could reconstruct f{uj) directly from the correlation function D{t) using data at all r, all the 
difficulties mentioned above would be avoided. Since the number of data for D{t) with a discrete set 
of time r is much smaller than the number of degree of freedom necessary for the reconstruction of 
/(w) in general, however, the standard fit is ill-posed for this problem. With some assumptions 
on the form of the spectral function the fit may work, but this is essentially equivalent to the 
multi-exponential or more complicated fit to the correlation function. 

In condensed matter physics, the reconstruction of the spectral function in quantum Monte Carlo 
simulations has been attempted with the maximum entropy method (MEM) It has been also 
successfully applied for image reconstruction in astrophysics. The most important assumption in 
MEM is that a probability for spectral functions can be assigned for given data of D{t). Then 
MEM can numerically reconstruct the most probable spectral function, using the Bayes's theorem in 
probability theory, without any strong constraints on its form. Recently, this method has been tested 
in lattice QCD |^,|| and first interesting results for the spectral function have been obtained |^-||]. 

In this paper, we employ MEM to reconstruct the spectral functions of pseudoscalar and vector 
mesons from the correlation functions previously calculated on lattices with the spatial size about 
3 fm at four different lattice spacings in quenched QCD . From the spectral functions we extract 
masses and decay constants for excited states as well as for the ground state. While they agree with 
results obtained from the exponential fits to correlation functions, errors for excited state masses 
from the spectral function are much smaller than those from the multi-exponential fit, so that we 
can estimate masses for excited states in the continuum limit with reasonable errors. We also find 
evidence that some excited states are composed of fermion doublers. 

This paper is organized as follows. In Sec. ^ we summarize our implementation of MEM and present 
results from tests using mock-up data generated from a realistic spectral function. Some details of 
the lattice QCD data and parameters used in our MEM analysis are given in Sec. III. In Sec. 



we present our results for the spectral function, which show excited state peaks as well as the ground 
state peak. From the position and the area of these peaks we extract masses and decay constants, and 
compare them with those obtained directly from correlation functions. The continuum extrapolation 
is made for these quantities. In Sec. 0, we argue that some peaks in the spectral functions correspond 



to the state containing two doublers of the Wilson quarks. Our conclusions are given in Sec VI. In 
the Appendix technical details of MEM are collected. 



II. MAXIMUM ENTROPY METHOD (MEM) 
A. Implementation 

The existence of a probability distribution for a spectral function is a key assumption in the max- 
imum entropy method. Using this assumption one can obtain the most probable spectral function 
for given lattice data D and all prior knowledge H , such as /(w) > 0, by maximizing the conditional 
probability 'P[F\DH\, where 'P[F\DH] is the probability of F with the condition that D and H are 
given. Here F stands for the spectral function /(w). Using the Bayes's theorem in probability theory 

Q. 

where P[^] is the probability of an event X, one rewrites the conditional probability V[F\DH] as, 

V[F\DH] cx V[D\FH] V[F\H]. (4) 

Here P[£'|FiJ] is the probability of data for a given spectral function, called the likelihood function, 
and P[i^|i7] is the probability of the spectral function for a given prior knowledge, called the prior 
probability. 
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The likelihood function is equivalent to in the least square method ||T^. For a large number 
of Monte Carlo measurements of a correlation function, the data is expected to obey a gaussian 
distribution according to the central limit theorem, which gives 



P[D\FH] = — 



Zl 

No 



(5) 
(6) 



with the normalization constant, = (27r)^°\/det C, and the number of temporal points No- The 
lattice propagator data averaged over gauge configurations, D{t), and the covariance matrix, C, are 
defined by 



^ * conA 



Cij — 



con! 

1 



NconfiNconf ^ 1) 



- J2 (Din) - D-in)) {D{r,) - D-ir,)) . 



(7) 



(8) 



where Nconf is the total number of gauge configurations and D'^{t) is the data for the n-th gauge 
configuration. Finally, Df{T) is the propagator constructed from the spectral function /(w) and the 
kernel K{u!,t) as 



Dfir) = / dujf{Lu)K{LU,T). 



(9) 



The prior probability is written in terms of the entropy S{f) [p^-p^ for a given model m(w) 
represented by a real and positive function, and a real and positive parameter a. The entropy S{f) 
becomes zero at its maximum point where /(cj) is equal to m{uj). Explicitly we have 



P[F\Hr 



1 



„aS 



S{f) 



Zs{a) 
duj 



E 

1=1 



f{uj) - to(w) - /(w)lo, 

fl 



/(^) 

m{u!) 



fi-mi- fi log 



mi 



(10) 
(11) 
(12) 



with the normalization constant, Zs{a) = (27r/a)^'^/^, calculated in Appendix [c|. In ( p2| ) the contin- 
uous spectral function f{uj) is approximately represented by a discrete set of points f{u>i) — fi with 
I = 1, ■ ■ ■ , N^. Hereafter we replace the prior knowledge H in by Hma, writing m and a explicitly. 
It is worth mentioning that this form of the entropy leads to a positive spectral function in MEM. 
Combining (|) and (|l^), one obtains 



P[F\DHma] cx 



ZLZs{a) 



QM) = aSif) 



(13) 



Therefore the condition satisfied by the most probable spectral function for a given a (and model 
171(10)) is given by 



Uf) 



Sf 



= 0. 



(14) 



The parameter a dictates the relative weight of the entropy S{f) and L. One can deal with a 
dependence of fa as follows. One first defines F[a\DHm] [^;^14|, the probability of a for given data 
and all prior knowledge, which can be transformed as 
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V\a\DHm]':xV[a\Hm] VF — . (15) 

J ZLZs[a) 

See Appendix for details. In the final result /{uj), a is averaged with this weight factor P[Q!|Di?m], 

f{uj)^ JdaP[a\DHm]faiuj)/ JdaP[a\DHm]. (16) 

This procedure is called Bryan's method |l^ and used in this article. We restrict the range 
of a in the actual average as amin < « < ctmax, where amin and amax are chosen to satisfy 
P[a\DHm] > lQP[amin.max\DHm] with a being the maximum value of P[a\DHm]. The standard 
choice of P[a|i?m] in ([l5| ) is either a constant or l/a [^I jl^jl^ . In the next section we will show that 
the final result is insensitive to the choice as long as P[Q;|I?-ffTO] is sharply peaked around a, and 
therefore we adopt P[Q;|i7m] — constant in our main analysis. 

In MEM it is not possible to assign error bars to each point in the spectral function since the errors 
between different points are strongly correlated. Instead we estimate the uncertainty of the spectral 
function averaged over w in a certain region by the method explained in Appendix I]. The magnitude 
of this uncertainty gives an estimate for the goodness of the given model m{uj) 



B. Test 

Several tests of MEM have already been carried out in Ref . |^ , where the dependence of results 
on the number of time slices N]j, the size of errors of data and the model m{uj) have been examined 
using mock-up data created from test spectral functions. The following conclusions have been drawn 
from the tests: 

1. Decreasing the error of data D{t) is more important than increasing No for obtaining better 
estimates of /(w) which reproduce the original spectral function more closely. 

2. It is better to include information of f{uj), such as the asymptotic value, if it is known, into the 
model m(Lu). 

3. If the obtained /(w) depends strongly on the model, a better model in the sense of leading to a 
/(cij) which is closer to the original spectral function gives smaller errors for the averaged f{uj). 

4. The error of the averaged f{ijj) in a certain region can be used to measure the significance of 
/(cj) in the region. For example, if the error of the averaged /(i^) around a peak is much smaller 
than the averaged value, the peak is likely to be true, and vice versa. 

Before applying MEM to actual data, we perform further tests on (l)thc dependence on No and 
the temporal separation of data Ar, and (2)the dependence on the choice of P[a|i?m]. For these 
tests we use a realistic spectral function in the vector channel of the e^e~ annihilation [|6|jT^], which 
is given by /m(w) = pin{oj)oj'^ , where the factor is expected from the dimension of meson spectral 
function, with 



TT 



'p2 rp(tj)mp 1 



(17) 



Here Fp is the residue of p meson resonance defined by 

{0\djf,u\p) = V2Fpmpef, = \/2/pm2e^, (18) 

with the polarization vector e^, and Tp{uj) includes the 0- function which represents the threshold of 
p ^ TTTT decay as 

1 ml / 4m2 \ ^ 
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We make dimensionful quantities dimensionless using the lattice spacing a, — > Lua, t r/a where 
a is set to 1 GeV^^. The values of parameters are 

mp = 0.77, = 0.14, Fp = 0.142, t^o = 1-3, S = 0.2, = 0.3, (20) 

where as is independent of uj for simplicity. The shape of Pin{i^) for this choice of parameters is shown 
in Fig. |l|. The value in the figure represents the area of Pm(w) for < a; < 6. 

We make mock-up data D{t) from /i„(a;) as follows, (i) The central value of D{t) is given by 
integrating the spectral function fini^j) and a kernel K(lu,t) = e~"^ over to in the same way as 
Df[T) in (|^). (ii) Errors of D[t) are generated by gaussian random numbers with the variance 
(T(Ti) = b ■ e°''^'D{Ti), a — 0.1, b — 10~^°, in order to incorporate the fact that the error of lattice 
correlation functions increases as r increases. 

In this test, no correlation between different r is taken into account, thus the covariance matrix 
C is set to be diagonal. The model function is given by m(w) = moo;^ with mo = 0.0277, which is 
motivated by the value of oo). We set the maximum value of w, uJmax = 6, and the w-space 

is discretized with an equal separation Auj — 0.01, and N,^ = 600. We also calculate the area of the 
MEM result Pout{^) for < a; < uJmax and define r = 'Y^i^i{pin{^i) ~ Pont(w/))^, to measure the 
difference between pin and Pout- 

We summarize the result of po«t(w) in various cases as follows. 
(1) To investigate the dependence of poMt(w) on Ar and N^^ we extract Pout{^) by MEM, from 



Data at large r are necessary to 
lat a false peak sometimes appears 



data with At = 0.5,0.33 and Njj — 16,31,46, as shown in Fig. 
reconstruct pout{^) at small oj correctly, as seen from the fact t 
around oj — from data with Ar = 0.5 and Ne, — 16 {rmax — At{Nd — 1) = 7.5) or with At = 0.33 
and N]j = 31 {rmax = 10). Once Tmax becomes large enough (larger than 15 in this case), a smaller 
At is better for the result, as seen from the comparison between results from data with At = 0.5 and 
At = 0.33 at No = 46. 

(2) We also check the dependence of Pout{'-^) on two forms of P[Q!|_ffTO], either P[a|iJm] constant or 
1/a. As shown in Fig. |[ the two choices give almost identical shapes of pouti^^), though the weight 
factor P[a\DHm] is rather different between the two cases. 

Our investigations add further information on the parameter dependence of the result in MEM, 
which we summarize as the following three points: 

5. T„iax — At{Nd ~ 1) must be sufficiently large for a reliable result of f{uj). 

6. Once Tmax is taken large enough, smaller At is better. 

7. The result Pout{^) is insensitive to the choice of P[Q;|i/TO]. 



III. LATTICE QCD DATA AND PARAMETERS IN MEM ANALYSIS 

We now apply MEM to the lattice correlation functions previously obtained in quenched QCD P,|l0| 
with the plaquette action for gluons and the Wilson action for quarks. The simulation was performed 
at four values of /3, corresponding to = 2-4 GeV for the continuum extrapolation, on 32'^ x 56 to 
64'^ X 112 lattices with the spatial size about 3 fm. Simulation parameters are compiled in Table |. 
At each /3, five values of the hopping parameter k, which correspond to m.,j /nip k, 0.75, 0.7, 0.6, 0.5 
and 0.4, were employed for the chiral extrapolation. The values of hopping parameters are numbered 
from heavy to light in Table |. For example, we call n corresponding to the lightest and the heaviest 
quark masses as i^51. Except for an additive renormalization factor, the average quark mass is equal 
to the average inverse hopping parameter given by 

K-'^\{^^' + ^2'). (21) 

where ki and K2 are the hopping parameters of quark and anti-quark in the meson. 

In our MEM analysis, we employ pseudoscalar and vector meson correlation functions, defined by 

^(Jr7/(T,x)(Jru)t(0,O)) = ( dujf{u:)K{uj,T), (22) 
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where T is 75(7^) for pseudoscalar(vector) meson, /(w) is a spectral function and K{uj, r) is a kernel. 
We use only point source data to satisfy the condition that f{oj) > 0. Since the spectral function of 
the meson propagator has dimension 2, we define a dimensionless function p{uj) as 

/H^pHa;^ (23) 

The model is chosen to be raiuj) = rriQUj'^ and the value of mo is taken equal to the asymptotic value 
of p{uj) in perturbation theory |6| given by 

/ n \ 

mo 



47r2 



where as is the strong coupling constant, the coefficients C^'s are perturbatively calculated in contin- 
uum QCD ||l9|, and Z is the renormalization constant for pseudoscalar (PS) or vector (V) operator. 
The spectral function from our data is insensitive to the value of too, as shown in Fig. y, where /(w) 
obtained with three different models are plotted for pseudoscalar and vector mesons at /3 = 6.47 and 
Kll. In the figure the horizontal bars indicate the region over which the result is averaged, while 
the vertical bars indicate the uncertainty in the averaged value of the result. Both averaged spectral 
functions and their uncertainty are almost identical for different models. Because of this property, we 
simply take as — 0.21 and employ the non-perturbative Zy and the perturbative Zps calculated at 



/3 = 5.90 in (24) for all /3. The normalization factor 1/2k is used also for the pseudoscalar meson with 
tadpole-improved Zps- Values of Z's as well as C^'s are given in Table ||. 

Other parameters in the MEM analysis such as Nd and {uja)max, are determined as follows. We 
take Nd as large as possible unless the error of data becomes too large for a reliable result, and we 
choose {Loa)max ^ and increase it until the result becomes stable. Both parameters are also given 
in Table For Aw, which should be smaller than 1/T, we take Aw = lO^"' around the peak of the 
ground state to determine the ground state mass accurately, while Aoj ~ 2.5 x 10^'^ away from the 
peak. 



IV. RESULTS 



In this section, we present our results of spectral function for pseudoscalar and vector meson prop- 
agators, from which we extract physical quantities such as masses and decay constants. 



A. Spectral functions 

Our results of p(w) obtained from meson propagators by MEM for three different K^^ at all j3 
are compiled in Fig. 0. The lowest peak corresponds to the ground state, the next peak corresponds 
to the first excited state and so on. At fixed /3, positions for these peaks move toward smaller u as 
the quark mass decreases. This shows that meson masses decrease with decreasing quark mass, as 
expected. The number of peaks increases from (3 — 5.90 io (3 — 6.47 for both pseudoscalar and vector 
channels, since more states with higher energy appear in spectral functions for larger lattice cut-off 
(smaller lattice spacing). All peak positions move to smaller values as (3 increases, except the peaks 
at uja « 1.7 for the pseudoscalar channel and at oja « 2 for the vector channel. Thus masses in the 
physical limits stay finite, except those of the latter peaks which become infinite. We discuss these 
unphysical states in more detail in the next section. 



B. Meson masses 



From peak positions of the spectral function, we determine masses of excited states as well as the 
ground state. Errors of these masses are estimated by the single elimination jackknife method. 

In order to check whether the peaks in the spectral function really correspond to particle states 
in correlation functions, we also extract masses of the ground and the first excited states by fitting 
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correlation functions with a double exponential form. In order to obtain the mass of the first excited 
state reliably, correlation functions from both point and smeared sources for the p meson are simulta- 
neously fitted. Results at /3 = 5.90 are given in Table III and Fig. ^, where errors are again evaluated 
with the single elimination jackknife method, together with those obtained by MEM. We find that 
the ground state masses from the two methods agree very well, and the first excited state masses are 
consistent with each other within the statistical error. It is noted, however, that the error for the first 
excited state by MEM is much smaller than the one by the usual fit. 

We determine the chiral limit and the critical hopping parameter Kc where the ground state of the tt 
meson mass vanishes by extrapolating (mTrCt)^ linearly in K~^. For other states including the excited 
states of TT mesons, the masses ma themselves obtained from the spectral function are extrapolated 
linearly in to the chiral limit. The chiral extrapolation at each P is shown in Fig. 0. Some 

excited state peaks do not appear in the spectral functions obtained from some jackknife samples. 
These masses are excluded from the chiral extrapolation and are not plotted in the figures. The 
lattice spacing a is fixed by setting the ground state mass for the p meson in the chiral limit to the 
experimental value, rrip = 770 MeV. All dimensionful quantities are normalized by the p meson mass 
in the chiral limit. 

Masses in the chiral limit are compiled in Table IV , together with the result of the standard anal- 
ysis [p|jlC|] for the lattice spacing, which agrees with the values from the present MEM analysis. At 
f3 — 6.47, our lattice spacing has a larger error. This is caused by large errors of point source data 
at this f3. As shown in Fig. |8| the ground state masses for each K^^ agree with the previous results 
from exponential fit of smeared source data . 

Masses of excited states in the chiral limit are extrapolated to the continuum limit, except unphysical 
states mentioned before, as shown in Fig. ^. We see that the mass of the first excited state is consistent 
with the one reported in Ref. Q for both tt and p mesons. Note that the error for the first excited 
state of p meson from the double exponential fit at /? = 5.90 (square) is too large for a reasonable 
continuum extrapolation. Mass ratios in the continuum limit are given in Table 0. The mass ratio 
of the first excited state to the ground state for the tt meson in the continuum limit is 0.86(77), 
which should be compared with the experimental value 1.68(12), while the ratio for the p meson is 
2.00(74) in comparison to the experimental value of 1.90(3) or 2.20(2) (there are two candidates for 
the first excited state of the p meson in experiment). The first excited state masses for both mesons 
are consistent with experimental values albeit the errors are quite large. For the p meson we are not 
able to decide whether the first excited state is p(1450) or p(1700) due to the large error of our result. 



C. Decay constants 



From the spectral function we can also extract the decay constant for the ground state of tt and p 
mesons, /^r and fp, defined by 



(0|(J75^.)^a*|7ro,p = 0) = ^^^I^II^±f[ ^ , (25) 

{0\{dlf.u)iat\po,P ^ 0) ^ V2fpm%—l[J — . (26) 



We employ the one-loop result with tadpole-improvement for the renormalization factor Za |2^ given 
by = 1 — 0.316ay(l/a), and the bare quark masses (m„ -|- md)^^^ derived from the axial ward 
identity |]lO|. For the vector meson decay constant, we use a non-perturbative value for JlOt . 

Decay constants can be extracted from the correlation function as follows. For the pseudoscalar 
meson we have, 

V(0|J75^^(T,x)(d75ii)t(0,0)|0) - ^(0|d75uK„)(^«|(rf75^^)^|0)^-^ (27) 

X n 

— . |(^oM75"|0)P^-^, r^cx), (28) 
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where En is n-th excited state energy, and a similar expression for the vector meson. Under the 
assumption that the ground state peak of spectral function is sharp, these correlation functions are 
related to the area of the spectral function around the ground state peak according to 

/ dco p,s{uj)uj^ = ^ 2^17 f ), (29) 

[ du;p.{oj)cu^ = f^ml^ f[ (^) . (30) 

J ground \ / 

For the first excited state, we also extract decay constants from the area of the spectral function 
around the first excited state under the same assumption as for the ground state. 

Decay constants obtained from the above relations are extrapolated linearly in K^^ to the chiral 



limit, as shown for P — 5.90 in Fig. |10|, and results are also given in Table [V. The decay constant for 
the first excited state of the tt meson should vanish in the chiral limit according to (^9|) , because 
remains finite in the chiral limit for excited states. This property can be seen from the figure. 

The continuum extrapolation is shown in Fig. ^l], and the results in the continuum limit are compiled 



in Table VI. For the ground state, the decay constant for tt and p mesons are consistent with previous 
results (squares) In the continuum limit we find /,ro — 80.3(5.9) MeV, which is smaller than 

the experimental value 93 MeV, and /p„ = 0.2062(84), which is slightly larger than the experimental 
value 0.198(4), and the first excited state decay constant for p meson /p^ = 0.085(36). 

D. Remark on spectral widths 

The width for the ground state peak should be zero for the tt meson, and should be very small for 
the p meson in the quenched approximation. Therefore the width for the ground state in spectral 
functions, if non-zero, is likely to be an artifact of MEM. The width F of the ground state peak for tt 



and p mesons are extrapolated to the chiral limit, and are compiled in Table IV. As shown in Fig. |1S_ 
these widths are very small and almost consistent with zero within errors, as expected. 

On the other hand, other states have larger widths. At this moment it is difficult to conclude 
whether these widths are physical or artifacts of MEM. In order to decide the nature of these widths, 
further researches are needed. 



V. UNPHYSICAL STATES AND FERMION DOUBLERS 

As mentioned in the previous section, the state in the pseudoscalar channel at uja « 1.7 and the one 
in the vector channel at ua « 2 appear with a large width in the spectral functions at all p. A similar 
state has been also observed in the Wilson quark action at /3 = 6.0 (a^^ — 2.2 GeV) of the plaquette 
gauge action |^ and at /3 = 4.1 {a~^ = 1.1 GeV) of a tree-level Symanzik improved gauge action 
We consider this state to be unphysical since its mass diverges toward t he co ntinuum limit. In fact 



the mass of this state can be fitted by Ci/a + C2 in Fig. h3 (see Table VII for numerical details), 
together with a linear continuum extrapolation for physical excited state. We also see from this figure 
that no physical excited states appear in the spectral function if its mass is larger than that of the 
unphysical state. At first sight, the state at uja w 1 seems to be a candidate of another unphysical 
state. We think, however, that this state is physical, since the position of the peak moves as /3 varies 
and moreover such a state has not been observed at a different lattice spacing . 

We argue that the unphysical state is a bound state of two fermion doublers of the Wilson quark 
action as follows. The pole mass of a free quark with Wilson parameter r = 1 is given by 

M{n) ^ -\og{l + ma + 2n) n = 0,1,2,3, (31) 
a 

where n = corresponds to the physical quark, and n 7^ represent doublers with n of the 3 spatial 
momenta components equal to 7r/a. At r = 1 the time doubler does not propagate due to its infinite 
mass. In the chiral limit the mass for the n = 1 doubler is given by M{l)a « 1.1, therefore, in this 
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free case, the mass of two n = 1 doublers is 2 x M(l) a « 2.2. Note that, for the meson correlation 
functions with zero spatial momentum, states consisting of, e.g., a physical quark and a doubler cannot 
contribute. 

In the interacting case, the mass for the bound state made of two doublers is expected to decrease 
from 2.2 in the free theory due to binding energy, which would depend on the quantum number of 
the state. This may explain the difference between the peak position at wa sa 1.7 for the pseudoscalar 
channel and at uja ~ 2 for the vector channel. 

From the consideration above we conclude that the unphysical state is a bound state of two n = 1 
doublers. We note that bound states of n > 2 doublers do not appear in the spectral function (in fact 
there is no peak at wa = 3.2 « 2 x M(2) a and 3.9 w 2 x M(3) a). 

VI. CONCLUSION 

In this study, we have applied the maximum entropy method to high-precision quenched lattice 
QCD data to extract the spectral functions for pseudoscalar and vector mesons. Masses for excited 
states as well as the ground state are obtained from the position of peaks in the spectral function, 
and decay constants are determined from the area under them. 

Masses of the ground and the first excited state agree with those obtained by the usual double 
exponential fit, showing the reliability of MEM, while the first excited state mass from the spectral 
function has much smaller errors, demonstrating the superiority of MEM. 

We have been able to make a continuum extrapolation for the first excited state for tt and p 
mesons, obtaining the masses 7717^ = 660(590) MeV and nip^ = 1540(570) MeV. While the errors are 
admittedly large, this is the first time that such an extrapolation has been attempted. For the ground 
state decay constant for tt and p mesons we have found that the result of MEM analysis is consistent 
with standard analysis. 

We have found a state in the meson spectral function at uja « 2 for all /3, and have argued that 
it is an unphysical bound state of two fermion doublers. If this interpretation is correct, this will 
be the first time that the doubler state has been identified numerically in lattice QCD simulations. 
Further confirmation of this interpretation can be made by changing the Wilson parameter r from 
unity, by analyzing the KS fermion data with MEM, or by considering meson correlation functions 
with a momentum of n/a. 

A future extension of MEM analysis is an application to unquenched data to see dynamical quark 
effects in the spectral function; decays and scatterings of intermediate states may be detected from 
possible widths in the spectral function. It will also be interesting to see the change of the spectral 
function before and after the phase transition at finite temperatures. 
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APPENDIX: TECHNICAL DETAILS OF MEM 

We collect technical details of probability theory and MEM in this Appendix. 

A. The Bayes's theorem 

In this section we list some results of probability theory and the Bayes's theorem used in MEM. 
The Bayes's theorem in probability theory flit] is given by. 
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where P[^] is the probabihty of an event X, and P[X|y] is the conditional probabihty of X given Y. 
These probabihties satisfy 

P[X]= J dYP[X\Y]P[Y], (A2) 
and the condition for normahzation, 

j dXP[X] = 1, (A3) 



J dXP[X\Y] = 1. (A4) 
In this article, we use P[X|yZ] which is the conditional probability of X given Y and Z. For P[X|yZ], 



(Al), (A2) and (A4) are rewritten, respectively, as. 



P[x|yzi . rasp, 

P[X\Z] = J dYP[X\YZ]P[Y\Z], (A6) 



j dXP[X\YZ] = 1. 



(A7) 



The most probable spectral function is obtained by maximizing the conditional probability P[F\DH] 
(in this section prior knowledge Hma is rewritten as H again for simplicity), and satisfies the condition, 

^M.o. (Aa, 

We rewrite P[F\DH] by the Bayes's theorem as, 

P[F\DH] = Pj^immH]^ 

The probability P \D \ FH] and P[F|i7] is the likelihood function and the prior probability, respectively. 



Integrating (|A9|) over F and using (A7), one finds that 



P[D\H] = J VFP[D\FH]P[F\H], (AlO) 

where VF is the measure of spectral functions. From this point of view, P[D|77] is a normalization 
factor related to the likelihood function and the prior probability, and we do not need to take account 
of it. 



B. Transformation of covariance matrix 



In this section we introduce the method which easily deals with a non-diagonal covariance matrix. 
If C is not a diagonal matrix, one can transform C into a diagonal from through C = Ra^R^^, where 
R is the transformation matrix and cr^ is the eigenvalue matrix of C . Kernel Ku = K{uji, Ti) and data 
Di = D{Ti) are transformed by R as. 
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Nd 

1=1 

i'=l 

After this transformation, the hkehhood function L defined in is written as, 

2 

J2 



i=i V 1=1 I 



This transformation does not require any changes in other parts of MEM. 



C. The normalization constant of the prior probability 



(Bl) 



(B2) 



(B3) 



The factor Zsipi) defined in (10) is the normahzation constant of the prior probabihty. In order to 
calculate Zsipi), we introduce a variable Xi which makes the curvature of S'(/) flat, and expand S'(/) 
by transforming /; into X; and applying the gaussian approximation to X{f) around X{m), 



S{f)^S{m)+Y,SXi 



1=1 



dS 
dXi 



X[m) 



- ^ 5Xi5X,, 



1,1=1 



dXidXi, 



X{m) 



U' ' ™ kk U 

where 6Xi = Xi{f) — Xi{m). From the property of Xi we choose 

dfi 



dX^, 



Since 



S{m) = 0, 



dS_ 



0, 



Ji 



we take the gaussian form for S{f), 



s{f)^--j2isx,r 



(Cl) 
(C2) 

(C3) 
(C4) 

(C5) 



(=1 



The measure VF is derived from the so-called 'Monkey Argument' [|6|, |13yiq ] and related to the metric 
of S{f). It is written as, 



dfi 



(C6) 



VF is transformed by (C3) such as, VF Y[i^=i dXi. We can easily integrate over // and obtain the 
normalization constant. 



Zs{a) = j VFe'^^^f^ 



lldXi 



exp 



1=1 



--aY^iSXi] 



(C7) 
(C8) 

(C9) 
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D. Uniqueness of MEM solution 



In this section we explain that the condition satisfied by the most probable spectral function, (A8) 
has only one solution, and has no local minimum. The likelihood function L satisfies 

H ^' af;jf =-X!3 -Q' with = (Dl) 
M'=i -1 ' '=1 

where zj's are non-zero real vectors and Zi^s are real vectors. The entropy and a real and positive 
parameter a satisfy 

where we have used < /; < oo and < a < oo. The matrix d'^Qa{f) / dfidfii is negative definite. 



Using RoUe's theorem, one can verify that ( |Aq ) has only one solution corresponding to the global 
maximum of Qa{f), if it exists dj. Roughly speaking, since the curvature of Qa{f) is always negative, 
Qaif) has only one maximum. 



E. The calculation of P[a\DHm] 



In order to search for the most probable value of a, we need to evaluate the conditional probability 
P[a\DHm]. This conditional probability is used in Bryan's method as the weight factor of 
averaging over a. In order to calculate P[a\DHm], we transform P[a\DHm] by the Bayes's theorem 
and rtA6|) as, 



P[a\DHm] = P[D\Hma]P[a\Hm]/P[D\Hm] 

= P [a I Hm] jvFP[D\ FHma] P [F \ Hma] /P[D\ Hm] 



CX P[q! 



\Hm] jvF- 



ZLZs{a) ' 



(El) 
(E2) 

(E3) 



Under the assumption that P[F\DHma\ is sharply peaked around the most probable spectral function 
fa, Qaif) is expanded in the variable Xi{f ) used in Appendix [c| and the gaussian approximation 
around Xi{f) = Xi{fa), 



Qa{f)~Qa{fa)+Y.^^l 



dQa 



1=1 



dX, 



Qa{fa)+J2^^' 



dfi- dQc 



n' 



dxi dfi. 



Ll' = l 



dXidXf ^(^^) 
- 2^ 5Xi5X, — -———^ 



kk'll' 



where 6Xi — Xi{f) — Xi{fa). Because 



dQa 



dfi 



0, 



d^Qa 



dfidfi' 



"a 
fi 



(E4) 
(E5) 



dfidfi'J 



(E6) 



we can write 
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i/=i 

where A;;' is a real symnietric x matrix defined as 



We then obtain 



^fl^f^ 



I' 



cx P[a|i7m] e'^"^^") J| 
here A; 's are the eigenvalues of A 



Ll' 



1 = 1 



a + Xi 



(E7) 
(E8) 

(E9) 
(ElO) 



F. Estimation of uncertainty in MEM 

In MEM, it is possible to estimate the uncertainty of spectral function averaged over a certain region 
I of 



(/«> 



(Fl) 



where {&) = / T>F QP[F\DHma]. Using the gaussian approximation and the variable Xi{f) in 
Appendix the covariance of the spectral function can be calculated as, 



(F2) 
(F3) 

(F4) 



where T — a5 + K. The form of (F4) is readily available because it is the Hessian of the Newton search 
algorithm P,|6|,[l7|| used to find fa- The uncertainty is estimated as. 



(F5) 



Similar to the spectral function, the error of averaged spectral function in a certain region I is averaged 
over a with the weight factor P[a\DHm], 



(Sfh 



J daF[a\DHm]^/{(SU? 
J daP[a\DHTn] 



(F6) 
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TABLE I. Simulation parameters of hadron propagator data |^|lO| used in the present MEM analysis. The 
numbering of hopping parameter is introduced for convenience. The smallest number corresponds to the 
heaviest quark mass, and vice versa. 



hopping parameter k 



/3 


lattice sizc{L-^T) conf.# 


swcep/conf. 


/3 


1 2 3 4 5 


5.90 


32^ 


56 


800 


200 


5.90 


0.1566 0.1574 0.1583 0.1589 0.1592 


6.10 


40^ 


70 


600 


400 


6.10 


0.1528 0.1534 0.1540 0.1544 0.1546 


6.25 


48^ 


84 


420 


1000 


6.25 


0.15075 0.15115 0.15165 0.15200 0.15220 


6.47 


64^ 


112 


150 


2000 


6.47 


0.14855 0.14885 0.14925 0.14945 0.14960 



m^/mp 0.75 0.7 0.6 0.5 0.4 



TABLE IL Parameters used in MEM analysis. The bottom table shows {Nd, {^i}a)max)- 





Ci 


C2 


Z 


PS 


3/2 


11/3 


0.728 


V 


1 


1 


0.536 





5.90 


6.10 


6.25 


6.47 


PS 


(20,4.0) 


(32,4.5) 


(32,4.5) 


(45,4.5) 


V 


(21,4.2) 


(30,4.8) 


(30,4.8) 


(30,4.8) 



TABLE III. Comparison of the MEM analysis with the double exponential fit, using the vector meson 
correlation function at /3 = 5.90. The symbol Kn\n2 expresses the quark mass used in the correlation 
function, ni and 712 being defined in Table |. 





exponential fit 


MEM 




ground excited x^/d.o./ 


ground excited 


Kll 


0.5093(11) 1.08(11) 0.220 


0.5094(16) 1.034(30) 


K22 


0.4784(12) 1.08(14) 0.359 


0.4789(20) 1.018(37) 


if 31 


0.4772(15) 1.08(14) 0.466 


0.4779(20) 1.020(36) 


Ki2 


0.4613(15) 1.07(15) 0.587 


0.4623(23) 1.009(40) 


fs:33 


0.4435(22) 1.03(19) 0.687 


0.4451(27) 0.997(44) 


KAl 


0.4668(23) 1.09(17) 0.638 


0.4678(23) 1.020(37) 


K42 


0.4505(22) 1.06(22) 0.750 


0.4519(27) 1.006(44) 


KU 


0.4214(43) 1.08(21) 0.890 


0.4218(43) 0.969(58) 


K51 


0.4622(20) 1.15(21) 0.771 


0.4630(25) 1.020(40) 


K52 


0.4460(32) 1.11(19) 0.872 


0.4469(30) 1.004(46) 


K55 


0.4107(37) 1.19(20) 1.191 


0.4080(65) 0.929(70) 
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TABLE IV. Results obtained from the MEM analysis at each (3. Lattice spacings from the standard 
analysis |9|Jlo| are also listed. 



p 


5.90 


6.10 


6.25 


6.47 


a[GcY-'] 


0.503(6) 


0.387(6) 


0.321(5) 


0.220(25) 




1.986(25) 


2.583(40) 


3.105(53) 


4.52(51) 


a~'[GeY] ^,10| 


1.934(16) 


2.540(22) 


3.071(34) 


3.961(79) 




0.159881(13) 


0.154985(12) 


0.152556(9) 


0.149809(7) 


TT meson 




2.02(31) 


1.30(44) 


1.82(62) 


1.40(45) 






2.61(51) 


2.79(23) 


3.95(64) 




4.00(25) 


5.86(38) 


6.84(29) 


10.6(1.2) 


/tto/^^Po 


0.1157(21) 


0.1148(26) 


0.1099(28) 


0.119(14) 




0.036(16) 


0.028(14) 


0.029(21) 


0.007(4) 


p meson 




2.46(19) 


2.63(47) 


2.48(32) 


1.59(67) 


mpjnip,, 




3.81(65) 


4.02(41) 


3.53(71) 


mpjmp^ 








6.3(1.0) 




4.69(14) 


6.79(21) 


7.76(30) 


11.7(1.3) 


fpo 


0.2037(20) 


0.2088(25) 


0.2015(32) 


0.178(34) 


I pi 


0.1133(46) 


0.076(34) 


0.102(15) 


0.120(40) 




0.032(19) 


0.014(7) 


0.008(5) 


0.022(15) 



TABLE V. Masses of excited states normalized by the ground state p meson mass for the tt and p mesons 
in the continuum limit. Available experimental values are also given. 





m^^ / mp„ 




mpjmp„ 


m,pjmp„ 


continuum limit 


0.86(77) 


5.4(1.6) 


2.00(74) 


3.2(1.8) 


x'/d.o.f 


0.514 


0.538 


0.726 


0.240 


experimental value 


1.68(12) 




1.90(3) or 2.20(2) 



TABLE VL Decay constants for tt and p mesons in the continuum limit and experimental values. 





/to 


fpo 


fpi 


continuum limit 


80.3(5.9) MeV 


0.2062(84) 


0.085(36) 


x'/d.o.f 


0.618 


2.18 


0.555 


experimental value 


93 MeV 


0.198(4) 





TABLE VIL Fit parameters and /d.o.f. of the unphysical state fit for tt and p mesons. 





TT meson 


p meson 




2.57(30) 


2.924(25) 




-1.05(78) 


-1.051(58) 


x'/d.o.f 


0.3158 


1.476 



16 



0.15 



0.10 



0.1543 



3 



0.05 



2 4 6 

CO 

FIG. 1. The input spectral function pin{ijj). The value in the figure is the area under the curve for 
< w < 6. 



At = 0.5 I 



At = 0.33 I 



Nd = 16 



0.1635 
r=3.729 



0.1541 
r= 0.02860 



A^r, =46 



0.15 






0.10 




0.1540 






r=0.04738 


0.05 


\ 




0.00 


i 





1/ 



01574 
r=0.5235 





01540 




r=0.03189 


ti 

i 1 










2 4 




CO 



01540 
r= 0.02278 



FIG. 2. The output spectral function potit(w) obtained by MEM for different At and Nd is shown by solid 
lines. The input pin{ijj) is shown by long dashed lines. The values in each figure are the area of Pout{'jj) and 
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0.15 



0.15 



0.10 



3, 



P [a\Hm] =constant 
r=0.02125 



0.05 - i 



0.00 




2 4 
CO 



i 

a 



0.10 - 



3, 



0.05 



0.00 



P[a\Hm]=l/a 
r=0.02046 



2 4 
CO 



V[a\Hrn\= constant 
P[aWm]=i/a 




FIG. 3. Influence of the choice of P[a|iym]. The left figure is for P[a|_ffm] = constant, and the right 
for P[a|iJm] = 1/a. The figure below shows the corresponding 'P[a\DHm] normalized to unity for which 
data with At = 0.33 and Nd = 46 is used. The input /9ira(a') is shown by the long dashed lines, and 
r = X^^" (pin(w() — pout{i^i)Y represents the difference from pin{oj)- 



PS 



2m„ 



0.1 m„ 



V 



2m„ 



0.1m„ 



ay^i ma (oa 

FIG. 4. Model(mo) dependence for pseudoscalar (PS) and vector (V) channels at /3 = 6.47 and if 11. 
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FIG. 5. Spectral functions at all obtained by MEM for different values of . On the left hand side the 
TT meson spectral function and on the right hand side the p meson spectral function are shown. The state at 
wo 2 is considered as unphysical since its position does not move with /3. 
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FIG. 6. Comparison of the p meson mass for the ground and the first excited state from the spectral 
function and the one from the double exponential fit. Circles are slightly shifted to larger . 
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FIG. 7. Masses and their chiral extrapolations at all On the left hand side the tt meson mass and on 
the right hand side the p meson mass axe shown. Circles, squares, diamonds and left triangles represent the 
ground, the first excited, the second excited and the third excited state mass, respectively. The state shown 
by up triangles is considered unphysical as discussed in the text. Open symbols stand for the values in the 
chiral limit. 
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FIG. 8. Ground state masses of the p meson obtained by different analyses at /3 = 6.47. Squares are 
slightly shifted to larger . 
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FIG. 9. Continuum extrapolation of masses of physical excited states. For the first excited state, open 
diamonds and triangles represent the experimental value, and that obtained by Asakawa et al. [5]. For the p 
meson the open square shows the result of the double exponential fit at /3 = 5.90. 
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FIG. 10. Chiral extrapolations of pseudoscalar and vector meson decay constants a.t f) = 5.90. 
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FIG. 11. Continuum extrapolations of pseudoscalar and vector meson decay constants. Open diamonds 
show experimental values for the ground state. Open squares represent the previous results from standard 
analysis ]1C[|. 
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FIG. 12. Widths for the ground state peak of tt and p mesons and their continuum extrapolation 
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FIG. 13. Combination for the excited mass fit and the unphysical state fit of tt and p mesons. 
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